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ABSTRACT 

A  xmified  constitutive  material  model  was  implemented  in  to  the  PAFEC  Finite 
Element  Package  Level  8.1.  The  modifications  allow  the  PAFEC  Package  to  perform 
time  dependent  plasticity.  This  includes  the  calculation  of  the  updated  material 
Jacobian  used  m  the  formulation  of  the  global  stiffness  matrix.  Automated  time 
stepping,  reduced  storage  requirements  and  structural  convergence  features  were  also 
included  in  the  code.  A  series  of  test  analyses  addressing  different  aspects  such  as 
element  t3rpes  and  loading  conditions  were  performed  and  shown  to  agree  with 
expected  results. 
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Executive  Summary 

The  conditions  under  which  the  RAAF  operate  the  F-111  aircraft  are  conducive  to 
cracking  in  the  structurally  critical  wing  pivot  fitting,  and  certification  testing  in  the 
Cold  Proof  Load  Test  (CPLT)  has  demonstrated  failture  in  the  same  region.  The  main 
area  of  concern  is  the  Fuel  Flow  Vent  Hole  Number  13  (FFVH#13),  in  the  wing  pivot 
fitting.  In  order  to  calculate  inspection  intervals,  a  detailed  knowledge  of  the  elastic 
and  residual  stress  fields  is  required  for  these  locations. 

Qassical  plasticity  solution  techniques  reveal  severe  limitations  in  representing 
material  behaviour  tmder  non-symmetric  cycHc  loading.  As  a  result  AMRL  researched 
and  developed  an  alternative  constitutive  material  model  which  would  provide  an 
improved  representation  of  the  residual  stress  field  after  a  CPLT  t3q)e  load  cycle.  This 
report  describes  the  modifications  made  to  the  PAFEC  Finite  Element  Package  in  order 
to  implement  the  unified  constitutive  material  model,  and  shows  tiiat  it  produces 
validated  results. 

The  solution  capability  provided  by  the  implementation  of  the  unified  constitutive 
material  model  into  PAFEC  will  be  used  extensively  in  the  F-111  Structural  Integrity 
Task  to  perform  the  detailed  plasticity  analyses  on  the  F-111  FFVH#13.  The  elastic  and 
residual  stress  fields  developed  wiU  then  be  used  as  an  important  ingredient  for  the 
determination  of  the  inspection  interval  for  this  critical  location. 
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1.  Introduction 

riassiral  techniques  of  modelling  material  creep  and  plasticity  have  acceptable 
accuracy  provided  the  inelastic  strains  are  kept  small.  However,  they  are  poor  in 
describing  material  behaviour  where  significant  inelastic  strains,  followed  by 
unloading  and  cyclic  re-loading  are  involved.  Several  structurally  significant  problems 
on  the  F-lllC  aircraft  involve  large  cyclic  plasticity  effects.  The  Cold  Proof  Load  TesF s 
(CPLT)  non-S3rmmetric  cyclic  loading  induces  large  residual  stresses  into  the  structure 
and  these  have  a  substantial  effect  on  fatigue  life.  This  accordingly  influences  the 
Royal  Australian  Air  Force's  (RAAF)  structural  maintenance  policies. 

To  overcome  the  limitations  of  conventional  material  models,  new  unified  constitutive 
models  have  over  the  years  been  developed  by  researchers.  In  1988  the  Aeronautical 
and  Maritime  Research  Laboratory  (AMRL)  examined  several  of  the  material  models 
described  in  a  NASA  report  [1].  The  model  chosen  to  tackle  large  cyclic  plasticity 
effects  was  developed  by  Prof.  D.C.  Stouffer  and  V.  S.  Battachar  at  the  University  of 
Cincinatti,  Ohio,  USA.  Original  work  at  AMRL  utilised  a  one-dimensional  FORTRAN 
program  for  this  model  developed  by  Prof.  D.C.  Stouffer  and  colleagues,  reviewed  in 
Reference  2. 

In  1990  the  first  version  of  the  PAFEC  code  implementation  of  a  unified  constitutive 
material  model  was  completed  [3].  Reference  3  describes  the  modifications  that  were 
required  to  be  made  to  the  PAFEC  Finite  Element  Program.  The  initial  implementation 
resulted  in  very  small  time  steps  and  required  large  disk  storage  areas  as  a  result  of  tiie 
numerical  integration  scheme  employed.  In  addition  the  global  stiffness  matrix  was 
not  updated  at  each  iteration  which  reduced  structural  convergence  rates.  Recently 
development  work  addressing  these  deficiencies  via  the  solution  of  the  constitutive 
differential  equations  and  the  material  Jacobian  determination  has  been  performed  by 
Trippit  and  Jones  [4,5].  In  1995,  AMRL  contracted  Monash  University  to  provide 
details  of  title  solution  technique  and  material  Jacobian  determination  and  effect  the 
installation  of  these  improvements  into  AMRL's  version  of  the  PAFEC  code  [9]. 

This  report  describes  the  modifications  made  to  the  PAFEC  Finite  Element  Program  to 
incorporate  the  new  unified  constitutive  material  model.  Updated  input  data 
definitions  (PAFEC  User  Modules)  are  described  in  detail  in  order  to  allow  users  to 
apply  the  constitutive  model  in  their  analyses. 
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2.  The  Unified  Constitutive  Model 


Unified  constitutive  models  have  advantages  over  classical  techniques  which  allow 
improved  representation  of  strain  rate  and  cyclic  response  effects  [6].  These  effects 
allow  the  calculation  of  time  dependant  variables  such  as  strain,  creep,  stress 
relaxation  and  cyclic  hardening/  softening. 


The  unified  constitutive  model  implemented  into  PAFEC  FE  was  developed  initially 
by  Ramaswamy,  Stouffer  and  Bodner  [1]  and  further  developed  by  Trippit  et  al  [4  and 
5]  and  Searl  and  Paul  [7].  It  consists  of  the  inelastic  flow  equation  and  the  rate  forms  of 
the  state  variables:  back  stress  and  drag  stress  Z.  In  this  formtdation  the  inelastic 
flow  equation  is  written  as: 


where  s'  is  the  inelastic  strain  rate,  Kj  is  the  second  invariant  of  the  deviatoric  over¬ 
stress  tensor,  Sy-Cl^,  and  D  and  n  are  temperature  and  strain  rate  dependent 
material  parameters.  Here  the  evolution  equation  for  the  back  stress  Qjj  is  defined  as: 

Djjj  -  f2Sjj  +  Qjj  (2) 

where  f2  is  a  material.  The  inelastic  back  stress  Qjj  is  calculated  by  integrating  the 
following. 


=  h  eeff-2% 
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where  6  ?  is  the  inelastic  back  stress  rate,  fin>ax  is  the  maximum  value  of  back  stress 

expected,  s'g  is  the  effective  inelastic  strain  rate  and  fi  is  a  hardness  related  material 
parameter. 


The  drag  stress  Z  is  given  by. 


Z  =  Zi  +  (Zo-Z|)exp(-ms^) 


(4) 


where  Zo  is  the  initial  value  of  drag  stress,  Zi  is  the  final  value  of  drag  stress  and  m 
defines  the  rate  at  which  Z  evolves  from  Zo  to  Zi.  The  effective  inelastic  strain  s^,  is 
given  by: 
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where  is  the  inelastic  strain  tensor. 
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Details  of  the  derivation  of  the  material  parameters  D,  n,  f|,  £2,  ^max'  ^0' 
D6ac  steel  can  be  found  in  Reference  7. 


3.  Existing  PAFEC  Plasticity  Algorithm 

To  solve  non-linear  plasticity  problems  the  current  PAFEC  plasticity  algorithm  uses 
successive  incremental  application  of  loads  with  an  iterative  approach  using  the 
Newton-Raphson  scheme  to  obtain  convergence.  The  material  is  assumed  to  behave 
elastically  before  yield,  according  to  Hooke's  law.  Yielding  is  deemed  to  have  occurred 
when  the  effective  stress  of  a  Gauss  point  traverses  the  boxmdary  of  the  yield  function. 

There  are  currently  three  5deld  criteria  available:  1)  the  Von  Mises  criterion,  2)  tiie 
Drucker-Prager  criterion  and  3)  a  user  defined  criterion.  PAFECs  constitutive 
equations  use  the  Prandtl-Reuss  relationships,  and  optional  hardening  rules  can  be 
selected  ranging  from  pure  isotropic  hardening  to  pure  kinematic  hardening.  The  user 
is  allowed  to  choose  between  the  Euler  or  Runge-Kutta-England  integration  schemes. 
At  the  onset  of  plastic  yielding  the  algorithm  performs  a  sub-iteration  process  in  order 
to  converge  to  die  best  estimate  of  the  stress,  and  then  proceeds  to  check  for  structural 
convergence  via  the  residual  forces  (several  different  convergence  criteria  are 
available).  If  structural  convergence  has  not  been  achieved  then  the  program  wUl 
iterate  again  with  modified  increments  in  strains.  These  modified  strains  are  estimates 
based  upon  the  balance  between  the  internal  and  external  work  and  are  calculated 
from  die  residual  force  vectors.  Once  convergence  is  achieved  the  program  moves  to 
the  next  load  increment. 

The  main  controlling  FORTRAN  routines  are  shown  in  Figiure  1.  A  brief  description  of 
the  function  of  each  routine  is  also  included  in  this  figure  and  the  detailed  algorithm 
was  previously  described  in  [3]. 


4.  Implementation 


4.1  Introduction 

The  major  difference  between  the  new  constitutive  equations  and  the  standard  PAFEC 
solution  is  the  inclusion  of  time  dependent  plasticity  which  can  describe  creep  effects. 
Other  features  such  as  automatic  time  step  load  calculation,  an  additional  numerical 
solver  for  the  constitutive  equations,  material  Jacobian  calculation  and  structural 
convergence  time  stepping  control  have  also  been  included. 

The  fundamental  change  required  to  PAFECs  original  algorithm  involves  the 
replacement  of  the  state  determination  routines  and  the  re-definition  of  the  load 
increment.  The  revised  FORTRAN  subroutine  calls  referred  to  in  this  report  are 
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presented  in  Table  1  (new  routines)  and  Table  2  (old  PAFEC  modified  routines).  An 
updated  program  flow  chart  of  the  modified  areas  has  been  provided  in  Figure  2. 

Specific  details  of  the  overall  implementation  may  be  found  in  Section  4  of  Reference  3 
and  are  only  be  repeated  here  where  changes  have  been  made.  A  description  of 
additional  and  enhanced  features  implemented  by  this  work  follows. 


4.2  Data  Input 

The  PAFEC  system  allows  users  to  define  data  modules  which  enable  easy  input  of 
data  into  tihe  program.  With  the  unified  constitutive  model  implementation,  four 
additional  data  modules  were  required.  The  new  data  modules  are: 

1.  AUTOMATIC.CONTROL 

2.  CONSTANTS.NEW.FLOW.LAW 

3.  UNIFIED.CONTROL 

4.  TOTAL.TIME.VERSUS.FORCE.DATA 

These  modules  and  their  influence  on  the  behaviour  of  the  solution  are  described  in 
Appendix  A.  The  PAFEC  modular  dictionary  has  also  been  modified  to  include  the 
above  modules.  Details  of  these  modifications  can  be  found  in  Appendix  B,  Section  B.l. 


4.3  Initialisation 

Modifications  made  to  the  initialisation  routines  are  detailed  below: 

4.3.1  Control  flags 

The  LFLAGS  common  block  has  been  extended  in  routine  R09702  to  include  positions 
30  to  44.  Details  of  the  meaning  of  each  new  control  flag  can  be  found  in  Appendix  B, 
Section  B.2. 

4.3.2  Data  validation  and  initialisation 

PAFEC  is  divided  into  10  Phases  which  perform  different  functions,  see  [8].  The 
PAFEC  system  generally  performs  all  data  validation  in  Phases  1  to  4,  but  the  new  data 
modules  are  validated  at  the  initialisation  stage  of  Phase  9.  Default  values  are  inserted 
where  applicable  if  none  have  been  provided,  and  WARNING  messages  are  printed  to 
the  output  files.  The  routine  JPIOOO  performs  this  function. 

4.3.3  Space  allocation 

In  the  previous  implementation  a  large  number  of  non-linear  data  sets  were  set  up  to 
handle  the  information  required  to  be  stored  on  a  Gauss  point  basis.  This  has  been 
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made  more  efficient.  Appendix  B,  Section  B.3  contains  details  of  the  changed  data  set 
configuration  with  information  for  the  NLSET  storage  array. 

4.3.4  Additional  scratch  modules 

Additional  scratch  modules  are  required  to  hold  solution  control  information,  double 
precision  material  data,  solution  recovery  information  and  other  data.  Details  of  these 
modules  can  be  foimd  m  Appendix  B,  Section  B.4. 


4.4  Multiple  Load  Case  Capability 

Code  in  the  previous  version  was  not  capable  of  handling  more  than  one  load  case. 
PAFEC  does  not  have  this  limitation  and  the  routines  were  re-written  allow  multiple 
load  cases.  The  input  of  total  time  data  via  the  TOTAL.TIME.VERSUS.FORCE.DATA 
module  (see  Appendix  A)  has  been  affected  by  these  modifications.  This  module  must 
now  contain  the  total  time  data  for  each  individual  load  case  corresponding  to  the 
INCREMENTAL  module.  New  routines  JP1105  and  JP1106  have  been  included  to 
process  the  Load-Time  data  and  the  old  routines  JPllOO  and  JPllOl  have  been 
modified.  A  negative  total  time  is  permitted  to  indicate  an  unload  condition  so  that  the 
elastic  stiffness  matrix  should  be  used  at  the  next  load  increment. 

The  capability  to  have  several  load  cases  active  at  any  one  time  has  been  incorporated 
into  the  modifications,  although  there  is  no  current  requirement.  The  routines  would 
require  validation  if  this  fadlity  is  needed  in  the  future. 


4.5  Load  Step  Housekeeping 

The  majority  of  the  work  required  to  be  performed  prior  to  each  load  step  has  not  been 
changed  in  this  version,  see  [3]  for  specific  details  of  previous  changes.  However, 
modifications  have  been  made  to  allow  for  the  overwriting  of  the  current  load 
increments  data.  This  allows  many  sub-increments  to  be  performed  while  storing  only 
selected  increments  at  specified  time  points.  This  feature  is  controlled  by  the 
DELTA.STORAGE.TIME.VALUE  option  in  the  ALTrOMATICCONTROL  module.  The 
BLOCK.SIZE  option  allows  the  program  to  expand  the  incremental  and  dynamic 
incremental  modules  beyond  the  original  size  defined,  thus  defining  the  complete  size 
of  the  stress/ strain  output  files.  Correct  use  of  these  two  options  can  limit  the  size  of 
disk  storage  requirements  for  a  particular  analysis. 

The  standard  increment  cormter  (IBASE(153))  has  been  replaced  with  a  sub-increment 
cotmter  (IBASE(180))  that  indicates  of  the  total  number  of  increments  a  solution  has 
performed.  The  standard  increment  cotmter  gives  the  number  of  increments  that  have 
actually  been  written  to  the  results  files. 
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4.6  Material  Jacobian  Determination  using  Perturbation  Theory 

Reference  9  shows  how  a  material  Jacobian  is  calculated  for  each  Gauss  point  widiin 
the  yielding  zone,  using  perturbation  theory.  The  UNIFIED.CONTROL  module  allows 
control  over  the  solver  method,  perturbation  parameters  and  solver  tolerances.  The 
routine  kmodel  is  the  main  routine  that  handles  the  unified  constitutive  solution  of  the 
equations. 

The  Jacobian  is  calculated  at  the  material  level  and  then  utilised  in  the  PL6XXX  series 
of  routines  when  the  tangent  stiffness  matrix  is  formulated.  Utilising  this  technique 
has  achieved  a  dramatic  improvement  in  the  number  of  iterations  required  to  achieve 
structural  convergence  within  an  analysis.  The  calculation  of  the  Jacobian  by 
perturbation  can  be  quite  expensive  in  terms  of  CPU  time,  however  this  is  more  than 
compensated  by  the  improved  structural  convergence  rates.  In  order  to  reduce  the 
penalty  of  calculating  the  Jacobian,  it  was  been  linked  to  the  FREQUENCY  and 
F1X.STIFFNESS  options  within  the  CONVERGENCE  module,  see  [8]. 

4.7  Stress  Determination 

Three  t3q)es  of  elements  are  currently  supported:  1)  3D  Isoparametric  elements;  2)  2D 
plane  stress  and  plane  strain  elements;  and  3)  Semi-Loof  shell  elements.  Each  main 
routine  has  been  modified  to  call  die  main  unified  constitutive  state  determination 
routine  JP3000  which  in  turn  calls  kmodel.  The  routine  kmodel  performs  the  Jacobian 
perturbation  and  kkmodel  performs  the  state  determination  for  the  unified  constitutive 
model.  In  the  previous  version  JP3000  called  JP4000  and  the  JP4XXX  series  of  routines . 
has  been  removed. 

Elements  that  are  defined  to  be  elastic  b}q>ass  the  call  to  kmodel  and  the  stress  is 
calculated  using  the  standard  techniques.  For  efficiency,  elements  defined  in  the 
YIELDING.ELEMENTS  modide  [8]  are  also  checked  to  see  if  plasticity  is  occurring  and 
if  not  are  returned  to  have  the  elastic  stress  calculated  using  the  standard  techniques. 

It  should  also  be  noted  that  the  strain  definition  varies  between  PAFEC  and  the 
constitutive  model  equations.  The  inelastic  constitutive  equations  use  the  tensor 
definition  of  strain,  but  results  are  output  as  engineering  strains.  PAFEC  works 
internally  wifii  engineering  strains,  but  outputs  its  results  in  terms  of  tensor  strain. 
This  has  been  taken  care  of  in  the  routine  JPSOOQ. 


4.8  Livermore  Solver 

The  stress  determination  routine  kkmodel  utilises  the  Livermore  solver  [12]  for  the 
ordinary  differential  equations  of  the  unified  constitutive  model.  This  solver  is  very 
good  when  used  with  stiff  and  non-stiff  systems  of  differential  equations.  It  has 
proved  to  be  very  successful  for  strain  rate  independent  materials  which  cause  the 
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equations  to  be  extremely  stiff.  The  previous  implementation,  discussed  in  [9],  had  to 
talce  very  small  material  time  steps  which  resulted  in  very  small  structural  time  steps 
and  large  solution  times  (in  the  order  of  weeks  real  time  for  the  typical  F-lllC 
problem).  The  Livermore  solver  performs  this  small  stepping  internally  and  very 
efficiently. 


4.9  Post-Convergence  Processing 

The  current  modifications  allow  for  the  same  convergence/non-convergence  handling 
outiined  in  [3].  The  only  modification  is  the  time  step  determination  technique  which 
has  been  re-worked  to  take  full  advantage  of  the  large  material  time  steps  tihat  are 
achieved  with  the  use  of  the  Livermore  Solver. 

4.9.1  Time  step  determination 

In  automatic  mode,  the  time  step  for  the  next  load  step  is  based  upon  the  ability  of  the 
material  law  to  maintain  integrity  as  well  as  how  the  structure  is  performing  in  terms 
of  global  convergence. 

The  unified  constitutive  routine  returns  a  delta  time  expansion  or  contraction  factor 
based  upon  the  level  of  plasticity  at  every  Gauss  point.  The  original  version's  routine, 
JP3920,  stiU  determines  the  minimum  time  step  for  the  structure  based  upon  this 
factor.  Previously  this  was  used  for  the  whole  structure,  however  with  the  improved 
solver,  the  time  step  tends  to  be  quite  large,  eg  it  is  not  unrealistic  to  be  able  to  jump 
from  load  turning  point  to  load  turning  point. 

The  routine  JPllOO  has  been  significantly  modified  to  compute  the  global  time  step 
based  upon  the  following  criteria: 

1.  material  convergence, 

2.  structural  convergence, 

3.  an  applied  load  turning  point, 

4.  turning  point  recovery  (time  step  after  a  turning  point  as  been  reached), 

5.  non-convergence  situations: 

a)  material  non-convergence, 

b)  maximum  iterations  exceeded. 

The  user  can  tailor  structural  convergence  parameters  in  the  AUTOMATIC. 
CONTROL  module  (see  Appendix  A  for  specific  details). 

4.9.2  Phase  Output  Control 

The  current  time  dependent  solution  algorithm  can  involve  many  time  steps,  and  as 
before  there  is  a  need  to  limit  unnecessary  I/O  in  all  areas  of  the  program.  Since  the 
last  report  [3],  modifications  have  been  made  to  improve  the  ability  to  extract  only  the 
required  information.  This  relies  heavily  on  the  ability  of  PAFEC  Integrated  Graphics 
System  to  extract  the  stress  and  strain  data  that  is  required.  Output  to  the  Phase  output 
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file  can  be  controlled  by  the  OUTPUT.PRINT.  CONTROL  option  in  the  AUTOMATIC. 
CONTROL  module  and  this  is  described  in  Appendix  A. 

The  table  output  option  has  not  been  changed  and  is  still  available. 

4.9.3  Reduced  Integration 

In  order  to  attain  faster  solution  times,  subroutine  NL0007  was  modified  so  that  all 
available  elements,  except  semi-Loof,  utilise  reduced  integration.  This  routine  has  been 
taken  out  of  the  global  library  and  now  has  to  be  used  separately  when  required. 
Using  this  option  can  significantly  reduce  the  size  of  memory  (BASE)  required  to  solve 
large  problems. 


4.10  Control  Option 

A  control  option  CONSITTUnVE  has  been  included  into  the  PAFEC  system  so  that 
the  unified  constitutive  law  can  be  utilised  at  any  time.  This  option  indicates  to  the 
PAFEC  system  that  the  unified  constitutive  object  code  is  to  be  linked  in  to  the  Phase  9 
executable. 

4.11  Convergence  Module 

The  CONVERGENCE  module  has  some  very  useful  features  that  can  improve  the 
performance  of  an  analysis.  Care  has  been  taken  to  make  sure  that  these  options  are 
active  and  working  correctly.  Note  that  the  load  increment  now  refers  to  the  stored 
load  increment  and  not  the  actual  number  of  increments  performed  (see  Section  4.5). 

4.12  Tolerance  Module 

The  TOLERANCE  module  hke  the  convergence  module  can  be  used  to  tailor  and 
improve  the  performance  of  an  analysis.  This  module  again  uses  the  stored  load 
increment  and  not  the  actual  number  of  increments  performed.  TOL5  has  been 
generally  used,  however  diis  tolerance  formulation  can  cause  difficulties  when  hying 
to  converge  on  a  zero  applied  load  value.  In  these  cases  the  TOL7  (not  documented  in 
the  PAFEC  Level  8.1  manual)  has  been  used  successfully. 

4.13  PAFEC  Modules  Disabled 

There  are  two  standard  PAFEC  plasticity  modules  ihat  are  only  required  so  that 
validation  may  be  passed,  but  have  actually  been  disabled.  They  are: 

1.  STATE.DETERMINATTON 

2.  UNIAXIAL.PROPERTIES 
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5.  System  Testing  and  Validation 

5.1  Introduction 

The  constitutive  material  model  was  validated  with  respect  to  the  PAFEC  system  to 
ensure  that  no  array  corruption  occurred,  stresses  and  strains  components  were  correct 
and  the  element  stiffness  matrix  calculation  was  correctly  updated  in  the  global 
stiffness  matrix.  Here,  the  results  generated  from  the  test  jobs  were  compared  directly 
with  an  external  FORTRAN  program's  results  for  the  same  set  of  material  parameters. 
The  validation  of  the  unified  constitutive  model,  with  respect  to  D6ac  steel  can  be 
foimd  in  References  7, 10  and  11,  where  numerical  results  were  directly  compared  to 
experimental  test  data. 


5.2  PAFEC  Options  Validated 
5.2.1  Description 

Two  and  tiuree  dimensional  test  cases  were  generated  to  demonstrate  that  the  material 
code  was  behaving  as  expected.  These  test  cases  included  all  the  element  types 
outlined  in  Section  5.2.2.  Combinations  of  multiple  elements,  element  types,  uniaxial 
and  cyclic  loading,  plane  stress  and  plane  strain,  and  additional  elastic  elements  were 
investigated  to  ensure  that  the  code  did  not  affect  or  cause  corruption  to  other  aspects 
of  the  finite  element  package.  The  generated  solutions  were  compared  against  known 
answers  for  the  particular  case  being  investigated.  Table  3  gives  a  summary  of  tiiese 
test  jobs  with  the  types  of  elements  and  loading  conditions  which  have  been  validated. 

This  section  presents  specific  results  for  two  of  the  validation  problems.  The  first 
problem  was  a  two  dimensional  plate  (Test  5,Table  3)  containing  4  plane  strain 
elements,  which  was  loaded  xmiaxiaUy  at  one  end  and  restrained  at  the  other,  see 
Figure  3a  and  3b.  The  structure  was  displaced  at  node  3  in  the  'y'  direction  and  all 
nodes  (119, 120, 121,  4)  in  the  same  plane  have  their  freedoms  repeated  to  provide  an 
even  distribution.  The  applied  loading  was  sufficient  to  induce  a  4%  strain  in  the 
material.  The  results  for  this  analysis  are  shown  in  Figure  4. 

The  second  validation  problem  was  a  three  dimension  bar  (Test  1,  Table  5)  which  is 
similar  to  the  previous  example  and  contained  8  isoparametric  brick  elements.  The 
displacement  loading  was  appHed  at  node  5  and  all  nodes  in  the  constant  'z'  plane 
have  their  freedoms  repeated,  see  Figmre  5a.  The  structure  was  restrained  on  three 
sides  as  shown  in  Figure  5b.  Again  the  loading  induced  4%  strain  and  the  results  are 
shown  in  Figure  6. 

The  simple  test  models  showed  that  mixing  the  element  types  produce  valid  results. 
Elements  with  a  local  axis  system,  different  from  the  global  axis  system,  have  been 
shown  to  produce  the  correct  global  stresses. 


9 


DSTO-TR-0529 


Simple  test  problems  cannot  always  trap  problems  such  as  array  overwriting.  A  more 
det^ed  validation  of  the  complete  implementation  of  the  unified  constitutive  model 
was  performed  for  a  circular  and  non-circular  hole  in  a  plate  [10,  11].  Here 
experimental  data  was  compared  to  the  finite  element  restilts  and  good  agreement  was 
achieved. 

5.2.2  Supported  Element  Types 

The  implementation  of  the  constitutive  material  model  supports  the  following 
elements: 

1.  3D  isoparametric  brick  and  prism  elements: 

a)  37100  (8  noded  brick) 

b)  37110  (20  noded  brick) 

c)  37200  (6  noded  prism) 

d)  37210  (15  noded  prism) 

2.  2D  isoparametric  rectangular  and  triangular  elements: 

a)  36200  (4  noded  rectangles) 

b)  36210  (8  noded  rectangles) 

c)  36100  (3  noded  triangles) 

d)  36110  (6  noded  triangles) 

Control  options  PLANE.STRESS  and  PLANE  STRAIN  have  been  tested. 

3.  2D  Semi-Loof  rectangular  and  triangular  elements 

a)  43210  (8  noded  rectangles) 

b)  43110  (6  noded  triangles) 

5.2.3  Convergence  Module 

Validation  of  the  CONVERGENCE  module  was  investigated  by  varying  parameters 
and  then  checking  that  the  program  performed  as  expected.  This  was  done  on  several 
of  the  test  jobs.  The  load  increment  defined  in  ihis  module  refers  to  the  stored  load 
increment  and  sometimes  it  is  possible  to  apply  tolerances  to  the  previous  increment. 
Care  must  be  taken  when  using  diis  feature. 

5.2.4  Tolerance  Module 

Validation  of  tire  TOLERANCE  module  like  the  convergence  module  was  performed 
on  various  test  jobs  and  foimd  to  work  in  conjunction  with  the  CONVERGENCE 
module. 
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6.  Conclusion 


This  report  has  outlined  the  implementation  of  the  unified  constitutive  model  into  the 
PAFEC  Finite  Element  Program  Level  8.1.  Several  new  enhancements  have  made  this 
implementation  capable  of  solving  large  cyclic  plasticity  problems  quickly  and 
efficiently.  This  version  has  reduced  solution  times  for  large  finite  element  models 
down  to  a  fraction  of  the  computational  time  that  was  previously  required. 
Application  of  this  work  will  enable  the  determination  of  the  residual  stress  fields  that 
exist  in  the  F-lllC  Wing  Pivot  Fitting  after  multiple  CPLT  load  cycles.  The  results  of 
that  analysis  will  provide  the  required  residual  stress  information  to  be  used  as  an 
essential  ingredient  in  the  determination  of  inspection  intervals  for  the  RAAF  F-lllC 
fleet  of  aircraft. 
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Tables 

Table  1:  Non  PAFEC  Subroutine  Description. 


Subroutine 

Name 

Description 

JPIOOO 

Initialisation  and  data  validation. 

JPlOlO 

Initialises  load  scratch  modules. 

JPlOll 

Extends  incremental  solution  size. 

JP1012 

Extracts  load  case  data  from  TOTAL.TIME  module. 

JPllOO 

Load  increment  control  routine.  Determines  structural  time  step. 

JPllOl 

Determines  %  load  based  on  current  time  for  each  Load  Case. 

JP1102 

Manipulates  incremental  and  djmamic  modules. 

JP1103 

Determines  if  current  increment  is  to  be  stored  to  the  results  files. 

JP1104 

Prints  the  summary  table  header  line. 

JP1105 

Interpolates  %  load  to  be  applied  time  for  each  load  case. 

JP1106 

Determines  if  an  elastic  solution  has  been  requested  at  a  turning 
point. 

JPlllO 

Determines  whether  material  Jacobian  is  to  be  calculated  for  use  next 
iteration. 

JP1120 

Reads/ writes  backup  of  non-linear  data  sets  and  displacements. 

JP3000 

Organises  variable  retrieval  and  storage,  executes  the  new  unified 
constitutive  algorithm.  Called  by  PL3000  and  PL3100. 

JP3002 

Extracts  unified  control  information. 

JP3100 

Internal  loads  for  thick  and  thin  shell  elements. 

JP3900 

Retrieves  stress  and  strain  components  for  a  single  Gauss  point. 

JP3901 

Retrieves  and  Stores  state  variables  for  a  single  Gauss  point. 

JP3910 

Stores  stress  and  strain  components  for  a  single  Gauss  point. 

JP3920 

Determines  minimum  material  time  step  for  next  load  step. 

JP3921 

Updates  reference  stress,  strain  and  state  variables. 

BT4800 

Calculates  equivalent  stress,  equivalent  plastic  strain,  equivalent 
elastic  strain  and  equivalent  total  strain. 

BT6000 

Controls  the  calctdation  of  the  elasto-plastic  modulus  matrix. 

BT6100 

Extracts  the  elasto-plastic  modulus  matrix  from  storage. 

JPUTIL 

Contains  some  commonly  used  routines. 

kmodel 

Main  unified  constitutive  routine.  Called  by  JP3000. 

kkmodel 

Performs  the  call  to  the  integration  scheme  for  each  time  step. 

kdy 

Provides  derivatives  as  a  function  of  time  for  Livermore  solver. 

kak2 

Calculates  K2  in  the  constitutive  equations. 

keffdvs 

Calculates  effective  stress. 

keffeps 

Calculates  effective  strain  rates. 

kepskk 

Calculates  the  sum  of  the  normal  components  of  the  total  strain. 

kwrkid 

Calculates  the  work. 

Isode 

Livermore  solver. 
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Table  2:  Modified  or  Re-written  PAFEC  Subroutines 


Subroutine 

Name 

Modifications 

R09702 

Sets  control  LFLAGS  in  common  block. 

NLIOOO 

Main  initialisation  control  routine. 

PL1020 

Retrieves  modules  228,229,230  and  231  into  BASE. 

NL1041 

Calculates  the  length  of  BASE  required  to  store  each  data  set. 

PL1040 

Creates  work  space  for  plasticity  solutions. 

NLllOO 

Main  load  increment  initialisation  control  routine. 

NLllOl 

Prints  the  load  summary  table  header. 

NLlllO 

Determines  the  type  of  solution  to  perform  next  iteration. 

PL3000 

Organisation  of  variable  retrieval  and  storage  prior  to  and  post  state 
determination  for  2D  and  3D  isoparametric  elements. 

PL3100 

Internal  loads  for  tiiin  and  thick  shell  elements.  Calls  JP3100. 

NL1150 

Writes  last  block  of  the  non-linear  data. 

NL1500 

Calculates  the  residual  forces. 

NL1600 

Tidies  up  at  the  end  of  an  incremental  solution. 

NL1300 

Checks  structural  sub-iteration  convergence. 

NL1340 

Determines  if  an  iterative  solution  has  converged. 

NL1360 

Print  and  nodal  extrapolation  control. 

R00411 

Stops  SS  file  warning  messages. 

R70620 

Stress  averaging  routine. 

NL0007 

Reduces  order  of  integration  to  2x2x2. 

L36300 

Controls  aU  stress/  strain  printed  output. 

L36301 

Controls  aU  stress/strain  output  to  results  files. 

PL5000 

Calls  relevant  yield  function  routine.  By-passed. 

PL6000 

Controls  the  calculation  of  the  elasto-plastic  modulus  matrix. 

PL6100 

Calculates  the  elasto-plastic  modulus  matrix. 

PL6200 

Calculates  the  elasto-plastic  modulus  matrix  for  2D  iso-parametric 
elements. 

NL9200 

Main  stress/strain  routine  for  361XX  and  362XX  series  of  elements. 

NL9220 

Calculates  and  prints  extrapolated  stress/strain  components. 

NL9222 

Calculates  principal  stresses  and  strains. 

NL9300 

Main  stress/ strain  routine  for  371XX  and  372XX  series  of  elements. 

NL9320 

Calculates  and  prints  extrapolated  stress/strain  components. 

NL9322 

Calculates  principal  stresses  and  strains. 

NL9400 

Main  routine  for  semi-Loof,  43XXX  series  of  elements. 

NL9410 

Sets  up  to  write  stress/strain  data  for  semi-Loof  elements. 

NL9412 

Writes  stress/strain  data  for  semi-Loof  elements. 

NL9500 

Main  stress/  strain  routine  for  beam,  34XXX  series  of  elements. 

NL9320 

Calculates  and  prints  extrapolated  stress/strain  components. 

NL9321 

Returns  stress/strain  values  at  comer  integration  points. 

NEWLIN 

Stops  new  lines  from  appearing  in  constitutive  output. 
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Table  3:  Summary  of  Test  Jobs^ 


Test 

# 

Data 

FUe 

Element 

Type 

Loading  Condition 

Description 

1 

T3D1 

37110 

Free  Uniaxial  Pull 

Test  the  general  stress  to  strain 
response  and  volumetric  response. 

2 

T3D2 

37110 

Constrained  PuU 

Test  the  general  stress  to  strain 
response  and  volumetric  response 
under  different  hydrostatic 

conditions  than  test  1. 

3 

T3D3 

37100 

Shear 

Test  the  general  stress  to  strain 
response  in  shear. 

4 

T3D4 

37110 

Free  Uniaxial  Cyclic 

Test  the  cyclic  stress  to  strain 
response. 

5 

TPSl 

36210 

Free  Uniaxial  Pull 

Plane  stress.  Test  the  general  stress 
to  strain  response  and  volumetric 
response. 

6 

TPS2 

36210 

Constrained  PuU 

Plane  stress.  Test  the  general  stress 
to  strain  response  and  volximetric 
response  under  different  hydrostatic 
conditions  than  test  5. 

7 

TPS3 

36200 

Shear 

Plane  stress.  Test  the  general  stress 
to  strain  response  in  shear. 

8 

TPS4 

36210 

Free  Uniaxial  Cyclic 

Plane  stress.  Test  the  cyclic  stress  to 
stratri  response. 

9 

TPEl 

36210 

Free  Uniaxial  Pull 

Plane  strain.  Test  the  general  stress 
to  strain  response  and  volumetric 
response. 

10 

TPE2 

36210 

Constrained  Pull 

Plane  strain.  Test  the  general  stress 
to  strain  response  and  volumetric 
response  under  different  hydrostatic 
conditions  than  test  9. 

11 

TPE3 

36200 

Shear 

Plane  strain.  Test  the  general  stress 
to  strain  response  in  shear. 

12 

TSHl 

43210 

Free  Uniaxial  Pull 

Test  the  general  stress  to  strain 
response  and  volumetric  response. 

13 

TMUl 

37110 

36210 

43210 

Free  Uniaxial  Pull. 

Mtdtiple  element  type  tests. 

^  Date  files  and  Result  files  reside  on  the  AMRL  computational  server  'bigted'  in  the  directory 
Vpafec/constitutive/source/validation/  runs'. 
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Create  workspace  modules  or  check  for  restart. 

Next  increment.  Re-init.  and  print  increment  summary. 
Next  iteration.  Check  for  updated  tangent  stiffness. 

Update  tangent  stiffness  and  internal  loads. 

Main  stiffness  calculating  routines  for  a  particular  Element. 
R37XXX  for  3D  Isoparmetric  Elements. 

R36XXX  for  2D  Plane  Stress  and  Plane  Strain  Elements. 
R43XXX  for  2D  Thin  Shell  Elements. 

Clean  up  after  stiffness  generation. 

Calculates  new  displacement  increments. 

Form  total  displacements  and  update  geometry. 

Prepare  for  residual  force  calculation. 

Sets  up  modules  prior  to  stressing. 

Main  stress  routines  for  a  particular  Element  T)^e. 

Cyclic  symmetry  module  swapping. 

Write  last  blocks  to  nonlioear  data  sets. 

Decide  which  method  to  use  to  generate  internal  forces. 
Update  tangent  stiffriess  and  internal  loads. 

Main  stiffness  calculating  routines  for  a  particular  Element. 
Qean  up  after  stiffness  generation. 

Update  any  element  loads  to  current  geometry. 

Update  tangent  stiffness  and  internal  loads. 

Main  stiffness  calculating  routines  for  a  particular  Element. 
Qean  up  after  stiffness  generation. 

Check  if  solution  method  is  iterative. 

Calculate  the  residual  forces  for  this  iteration. 

Print  solution  tolerance  metrics. 

Switch  the  ES  and  FS  files  if  necessary. 

Increment  has  converged,  or  iteration  Hmit  exceeded,  print. 
End  of  solution  house  keeping. 

Stress  averaging  performed. 

Write  last  block  of  stress  file. 


Figure  1:  Main  FORTRAN  plasticity  routine  calls 
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Initialisation 

R09702* 

NLIOOO*  -  PL1020* 

.  JPIOOO  -  JP1012 

-  JPlOlO  -  JP1012 
-  JPlOll 

-  NL1040  -  NL1041* 

-  PL1040* 


Load  step  house  keeping 

NLllOO*  -  JPllOO  -  JP1103 

-  JPllOl  -  JP1105 

-  JP1102 

-  JP1106 

-  JP1120 

-  NTLllOl* 

-  JP1104 

NLlllO*  -  JPlllO 

Stiffness  matrix  generation 

PL6000*  -  BT6000  -  BT6100 

PL6100*  -  BT6100 

PL6200*  -  BT6100 


Stress  determination 

PL3000*  -  JP3000  -  JP3900 

-  BT4800 

-  JP3002 

-  JP3901 

-  kmodel  -  kkmodel  -  Isode 

-  kepskk 

-  kepskk 

-  keffdvs 

-  kkmodel 

-  JP3901 

-  JP3910 

PL3100*  -  JP3100  -  JP3000  -  as  above 

Post  convergence  processing 
NL1150* 

NL1500* 

NL1300*  -  NL1340* 

NL1360*  -  JP3920  -  JP3921 

NL9XXX  series  of  stressing  routines 

Miscellaneous 

R00411*  L36300*  PL5000*  .  R70620* 

NL0007*  L36301*  BT4800  NEWLIN* 


NOTES 

*  =  Original  PAFEC  subroutine  modified. 

-  =  Subroutine  calls  this  routine. 

All  subroutines  beginning  with  JP  or  BT  have  been  written  from  scratch. 

Figure  2:  Program  Flow  Chart  of  Modified  Routines 
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Figure  4:  Results  for  Validation  Test  5. 
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Appendix  A 

INPUT  DATA  MODULE  DESCRIPTION 
A.1:  AUTOMATIC.CONTROL 

1  AUTO.CONTROL  |  TOTAL.T1ME.STOP  |  OUTPUT.FRINT.CONTROL  |  TABLE.PRIOT.CONTROL  | 
\  NUMBER.OF.LOAD.CYCLES  |  BLOCK.SIZE  |  DELTASTORAGE.TIME.VALUE  | 

\  REFACTOR  I  EX.FACTOR  |  RE.ITERATION  |  EX.ITERA1ION  | 

I  INITIAL.PELTA.TIME  |  ABORT.REPUCTION.TIME^FACTOR  |  MAXIMUMTIME.STEP  | 


This  module  defines  the  global  control  parameters  required  for  the  imified  constitutive  solution. 

AUTO.CONTROL  Specifies  whether  automatic  time  stepping  is  to  be  used. 

=  1  Automatic  control 

values  in  TOTAL.TIME.VERSUS.FORCE.DATA 

correspond  to  total  time. 

=  0  Manual  Control 

values  in  TOTAL.TIME.VERSUS.FORCE.DATA 

correspond  to  delta  time. 

THIS  OPTION  IS  RARELY  USED  AND  IS 
UNTESTED  IN  THIS  VERSION. 

Range:  0  and  1 
Default  =  2 

TOTAL.TIME.STOP  Determines  the  termination  time  for  the  solution. 

If  greater  than  the  maximum  referenced  time  in  the  TOTAL. 
TIME.VERSUS.FORCE.DATA  module  the  lower  value  is 
used  and  a  warning  message  is  printed  in  the  phase  output. 

Range:  >  0 
No  Default 

OUTPUT.PRINT.CONTROL  This  options  allows  the  output  written  to  the  phase  output 

file  to  be  controlled  or  totally  stopped. 


Value 

Meanine 

0 

No  output 

1 

Iteration  convergence  information  printed. 

10 

Displacements  printed. 

100 

Gauss  point  information  printed. 

1000 

Stress/ strain  data  printed 

10000 

Cumulative  reactions  printed 

2 

All  of  the  above. 
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Combinations  of  the  above  values  are  allowed  eg  if 
displacements  and  stresses  were  required  to  be  printed  then 
this  option  would  be  1010  ie  10  + 1000. 

Default  =  0 

TABLE.PRINT.CONTROL  Allows  a  table  of  loading  and  time  stepping  information  to 

be  written  to  the  phase  output. 

=  0  No  table  write 
=  1  Table  write 

Range:  0  and  1 
Default  =  1 

NUMBER.OF.LOAD.CYCLES  Specifies  the  number  of  times  the  loading  history  is  to  be 

repeated.  The  loading  is  defined  as  one  copy  of  the  force 
time  curve,  duplicated  for  each  load  case,  input  via  the 
INCREMENTAL  and  TOTAL.TIME.VERSUS.FORCE. 
DATA  modules. 

UNTESTED  IN  THIS  VERSION. 

Range:  ^1 
Default  =  1 

BLOCK.SIZE  This  is  the  total  number  of  increments  that  may  be  stored  in 

a  solution  run.  For  large  problems  this  may  cause  large  disk 
space  usage.  One  block  is  required  per  increment  stored. 

Range:  >  Maximum  number  of  increments  defined  in  the 

INCREMENTAL  module. 

Default  =  20 


DELTA.STORAGE.TIME.  VALUE 

This  is  the  time  interval  for  storage  of  results  in  output  files. 
Results  are  also  written  at  the  end  of  each  user  defined 
loading  point  irrespective  of  the  specified  value.  Results  can 
be  stored  every  increment  by  specifying  zero. 

Range:  ^  0 

Default  =  0 

RE.FACTOR  Specifies  the  factor  by  which  the  time  step  should  be 

reduced  if  the  iteration  count  from  the  previous  step  is 
greater  or  equal  to  RE.ITERATION. 

Range:  0  <  RE.FACTOR  <1 
Default  =  0.85 
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EX.FACTOR 


REITERATION 


EX.ITERATION 


INmAL.DELTA.TIME 


Specifies  the  factor  by  which  the  time  step  should  be 
expanded  if  the  iteration  coimt  from  the  previous  step  is  less 
or  equal  to  EX.ITERATION. 

Range:  ^  1 

Default  =  1.15 

Specifies  the  iteration  number  at  which  the  time  step  for  the 
next  increment  will  be  reduced  by  the  factor  RE.FACTOR 
See  Figure  A.l 

Range:  >  EX.ITERATIONS 
Default  =  8 

Specifies  the  iteration  number  at  which  the  time  step  for  the 
next  increment  will  be  expanded  by  the  factor  RE.FACTOR. 
See  Figure  A.1 

Range:  0  <  EX.ITERATIONS  <  RE.ITERATIONS 
Default  =  4 

Specifies  the  initial  time  step  at  the  start  of  nm  when  in 
automatic  control. 

Range:  >  0 
Default  =  2  seconds 


ABORT.REDUCTION.TIMEFACrOR 

Specifies  the  factor  by  which  the  time  step  should  be 
reduced  and  the  increment  re-attempted  in  the  case  of  a 
structural  convergence  failure,  maximum  iteration 
exceedance  or  a  material  point  convergence  failure. 


Range:  >  0 
Default  =  0.5 

MAXIMUM.TIME.STEP  Specifies  the  maximum  time  step  allowed.  This  is  used  to 

limit  the  time  step  from  expanding  past  unrealistic  values. 


Range:  >  0 
Default  =  5  seconds 
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Expand  dt 

Constant  dt  Reduce  dt 

1 

>  Abort 

1 

Expansion 

1 

Reduction 

Maximum 

Iteration 

Iteration 

Iterations 

Number 

Number 

Line  represents  number  of  iterations  to  convergence  for  the  previous  increment, 
dt  =  Time  Step 


Figure  A.l:  Structural  Time  step  calculation  logic. 
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A.2:  CONSTANTS.NEW.FLOW.LAW 


Defines  the  state  variables  properties  to  be  used  by  the  unified  constitutive  code. 

MATERIAL.NUMBER  Property  number  referenced  by  the  PLASTIC.MATERIAL 

module,  given  to  the  following  properties.  If  the  properties 
are  temperature  dependent,  the  MATERIAL.NUMBER  is 
repeated  for  each  temperature  where  the  state  variables  are 
known. 

No  Default 

TYPE.MATERIAL  Two  forms  of  the  unified  constitutive  equations  currently 

exist. 

=  1  General  purpose  equations 

=  2  Form  used  for  high  temperature  material  Hastelloy  X. 

Range:  1  and  2 
No  Default 

TEMP  Temperature  of  current  set  of  constitutive  state  variables 

properties. 

No  Default 

E  Young's  modulus. 

No  Default 

NU  Poisson's  Ratio. 

No  Default 

D  Inelastic  flow  scale  factor. 

No  Default 

N  Inelastic  flow  strain  rate  sensitivity  parameter. 

No  Default 

ZO  Initial  drag  stress. 

No  Default 

Z1  Saturated  drag  stress. 

No  Default 
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M 


FI 


F2 


OMMAX 


Cyclic  hardening  parameter. 
No  Default 

Strain  hardening  parameter. 
No  Default 
Yield  parameter. 

No  Default 

Maximum  back  stress. 

No  Default 
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A3:  TOTAL.TIME.VERSUS,FORCE.DATA 

I  LOAD.CASE  I  LIST  | 

Time  data  associated  with  the  INCREMENTAL  module  resulting  in  a  load-time  curve  to  be 
applied  to  the  structure. 

LOAD.CASE  The  number  of  the  load  case  in  which  the  TOTAL.TIME. 

VERSUS.FORACE.DATA  is  to  be  applied. 

7 

Default  =  lorl  +  the  previous  entry 

list  This  is  a  list  of  time  points  that  are  associated  with  the 

INCREMENTAL  module.  In  automatic  time  stepping  mode 
this  list  represents  the  total  time  to  obtain  the  required  load. 
In  manual  mode  this  is  a  list  of  delta  time  for  each  increment 
in  load  step. 

No  Default 

NOTE  1:  For  multiple  load  cases  a  value  of  0  is  acceptable. 

NOTE  2:  A  negative  value  infers  the  use  of  the  elastic 
modulus  at  the  next  increment. 

EXAMPLE  1  Load  Case 

TOTAL  TIME  VERSUS  FORCE  DATA 
LIST 

C - AUTOMATIC  MODE  — -  Total  =  30  seconds. 

10,  20,  30 

C  - MANUAL  MODE  —  Total  =2.4  seconds. 

.5,  .5,  .3,  .2,  .1,  .1,  .05,  .05,  .1,  .2,  .3 

EXAMPLE  2  Load  Cases 

TOTAL  TIME  VERSUS  FORCE  DATA 
LOAD  LIST 


C - AUTOMATIC  MODE  --  Total  =  60  seconds. 

1  10,  20,  30 

2  0,  0,  0,  40  50  60 
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A.4:  UNIFIED.CONTROL 


JACOBIAN.CALCULATE  |  METHOD.SOLUTION  |  ORDER.SOLUTION  | 
PERTURBATION.DISTANCEFACTOR  |  REEFSI  |  REOMI  |  REWRKI  j  REQEPsF 


AEEPSI 


AEWRKI 


AEQEPSI 


OOH.YOUR.STRAIN.IS.NOT.LINEAR  |  ELASTICJACOBIAN.FACTOR 


This  module  specifies  local  control  parameters  for  the  constitutive  model  including  solution 
schemes,  integration  tolerances  and  factors  controlling  material  Jacobian  determination. 


JACOBIAN.CALCULATE  Specifies  the  method  by  which  the  material  Jacobian  will  be 

calculated. 

0  =  Assumed  elastic. 

2  =  Numerically  determined  by  divided  differences. 
Perturbation 


Range:  0  and  2 
Default  =  0 

METHOD.SOLUTION  Specifies  the  scheme  for  the  constitutive  integration.  Refer  to 

Hindmarsh  [12]  for  full  description. 

=  10  Non  stiff  solution. 

=  21  Stiff  solution  requiring  full  Jacobian  matrix  df/dy. 

=  22  Backward  Differentiation  Formula  (BDF)  using 
numerical  Newton  requiring  full  Jacobian  matrix 
df/dy. 

=  24  Stiff  solution  reqtziring  banded  Jacobian  matrix  df/  dy . 
=  25  Backward  Differentiation  Formula  (BDF)  using 

numerical  Newton  requiring  banded  Jacobian  matrix 
df/dy. 

Range:  10,  21,  22,  24  and  25 
Default  =  22 

WARNING:  This  option  is  intended  for  experienced 
users  only. 

ORDER.SOLUTION  Specifies  the  maximum  order  for  the  constitutive  integration. 

Refer  to  Hindmarsh  [12]  for  full  description. 

Range:  Range  1-5  for  BDF  method  and  1-12  for  Adam's  method. 
Default  =  5 

WARNING:  This  option  is  intended  for  experienced 
users  only. 

PERTURBATION.DISTANCE.FACTOR 

This  field  specifies  the  perturbation  factor  which  controls  the 
perturbation  distance  (strain)  for  numerically  determined 
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material  Jacobians.  The  perturbation  distance  is  equal  to  the 
perturbation  factor  multiplied  by  an  estimate  of  the  error  on 
the  strain.  The  perturbation  factor  approximately  controls 
the  number  of  significant  figures  in  the  Jacobian.  For 
example  Fpert=1000  or  10^  will  result  in  approximately  3 
significant  figures  in  the  Jacobian.  Refer  to  [9]  of  this  report 
for  a  further  description. 

Range:  >  0 
Default  - 1000 


REEPSI 


REOMI 


REWRKI 


REQEPSI 


AEEPSI 


AEOMI 


NOTE:  This  control  is  coupled  with  integration  tolerances. 
An  increase  in  perturbation  factor  will  generally 
require  a  similar  decrease  in  integration  tolerances. 

Specifies  the  allowable  local  relative  error  on  the  integration 
of  inelastic  strain. 

Range:  >  0 
Default  =  5E-8 

Specifies  the  allowable  local  relative  error  on  the  integration 
of  inelastic  back  stress. 

Range:  >  0 
Default  =  5E-8 

Specifies  the  allowable  local  relative  error  on  the  integration 
of  inelastic  work. 

Range:  >  0 
Default  =  5E-8 

Specifies  the  allowable  local  relative  error  on  the  integration 
of  equivalent  inelastic  strain. 

Range:  >  0 
Default  =  5E-8 

Specifies  the  allowable  local  absolute  error  on  the  integration 
of  inelastic  strain. 

Range:  >  0 

Default  =  1E-11 

Specifies  the  allowable  local  absolute  error  on  the  integration 
of  inelastic  back  stress. 

Range:  >  0 

Default  =  2E-6 
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AEWRKI  Specifies  the  allowable  local  absolute  error  on  the  integration 

of  inelastic  work. 

Range:  >  0 
Default  =  lE-12 

AEQEPSI  Specifies  the  allowable  local  absolute  error  on  the  integration 

of  equivalent  inelastic  strain. 

Range:  >  0 
Default  =  lE-11 


OOH.YOUR.STRAIN.IS.NOT.LINEAR 

Specifies  an  allowable  nonlinear  strain  history  factor.  If  the 
strain  history  becomes  more  nonlinear  than  the  specified 
factor  a  warning  message  is  printed  to  flie  screen  and  to  the 
phase  9  output.  The  warning  messages  are  printed  when  the 
maximum  change  in  strain  rate  of  this  step  normalised  by 
the  maximum  strain  rate  of  the  previous  step  is  greater  than 
the  specified  factor.  For  example,  when  the  factor  is  set  to  0.4 
warnings  will  be  issued  when  the  change  in  strain  rate  is 
greater  than  40%.  This  option  is  useful  for  assessing  the 
accuracy  of  time  discretization  as  the  constitutive 
implementation  assumes  the  strain  history  is  linear  between 
solution  points.  If  warning  messages  occur  the  solution 
should  be  rerun  with  decreased  step  sizes.  This  option  can 
be  disabled  by  setting  the  factor  to  zero. 

Range:  ^  0 
Default  -  0  (disabled) 

Recommended  0.4  to  0.6 

NOTE:  This  field  controls  only  the  output  of  warning 
messages  and  does  not  affect  the  solution  in  any 
other  way. 

ELASTICJACOBIAN.FACTOR  Specifies  a  measure  of  inelasticity  below  which  the  Jacobian 

is  assumed  to  be  equal  to  the  elastic  stiffness  matrix. 
Numerical  perturbation  to  calculate  the  material  Jacobian  is 
only  performed  when  the  inelastic  strain  increment 
normalised  by  the  total  strain  increment  (at  the  material 
point  concerned)  is  greater  than  the  specified  factor.  The 
correct  use  of  this  factor  greatly  increases  the  performance  of 
near  elastic  material  points. 

Range:  >  0 
Default  =  0.05 
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Appendix  B 


B.l:  PAFEC  Data  Module  Definitions  For  The  Module  Dictionary 


CONSTANTS.NEW.FLOW.LAW 

MODULE  228 

Defaiilts 

MATERIAL.NUMBER 

-11 

1 

TYPE.MATERIAL 

-11 

1 

TEMP 

-11 

0 

E 

0 

0 

NU 

0 

0 

D 

0 

0 

N 

0 

0 

ZO 

0 

0 

Z1 

0 

0 

M 

0 

0 

FI 

0 

0 

F2 

0 

0 

OMMAX 

0 

0 

AUTOMATIC.CONTROL 

MODULE  229 

Defaults 

AUTO.CONTROL 

1 

1 

TOTAL.TIME.STOP 

0 

0 

OUTPUT.PRINT.CONTROL 

0 

1 

TABLE.PRINT.CONTROL 

1 

1 

NUMBER.OF.LOAD.CYCLES 

1 

1 

BLOCK.SIZE 

20 

1 

DELTA.STORAGE.TIME.VALUE 

0 

0 

RE.FACTOR 

0 

0 

EX.FACTOR 

0 

0 

RE.ITERATION 

8 

1 

EX.ITERATION 

4 

1 

INITIAL.DELTA.TIME 

0 

0 

ABORT.REDUCnON.TIME.FACrOR 

0 

0 

MAXIMUM.TIME.STEP 

0 

0 

TOTAL.TIMEVERSUS.FORCE.DATA 

MODULE  230 

Defaults 

LOAD.CASE 

-81  1 

LIST 

-40  0 

UNIHED.CONTROL 

MODULE  231 

Defaults 

JACOBIAN.CALCULATE 

0 

1 

METHOD.SOLUTION 

22 

1 

ORDER.SOLUTION 

5 

1 

PERTURBATION.DISTANCE.FACTOR 

0 

0 

REEPSI 

0 

0 

REOMI 

0 

0 

REWRKI 

0 

0 

REQEPSI 

0 

0 

AEEPSI 

0 

0 

AEOMI 

0 

0 

AEWRKI 

0 

0 

AEQEPSI 

0 

0 

OOH.YOUR.STRAIN.IS.NOT.LINEAR 

0 

0 

ELASTICJACOBIAN.FAC:rOR 

0 

0 
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B.2:  Control  Flags  Common  Block  Definitions 

Subroutine  R09702  sets  up  the  common  block  LFLAGS,  which  contains  a  number  of 
logical  flags  which  describe  the  current  analysis.  The  following  is  a  detailed  list  of  the 
common  block  LFLAGS  revised  from  Reference  3. 


NUMBER 

DESCRIPTION 

1 

2 

3 

4 

5 

6 

7 

8 

9 

10 

11 

12 

13 

INCREMENTAL  SOLUTION  FLAG. 

GEOMETRIC  NON-UNEARTTY  FLAG. 

LARGE  DISPLACEMENT  ANALYSIS. 

UNEAR  BUCKUNG. 

FREQUENCIES  OF  STRESSED  STRUCTURES. 

PLASnCTTY  ANALYSIS. 

CREEP  ANALYSIS. 

MATERIAL  NON-LINEARITY. 

SEISMIC  ANALYSIS. 

STRESS  ELEMENTS  IN  FRONT  ORDER. 

STRESSES  AT  INTEGRATION  POINTS. 

INCREMENTAL  STRESS/STRAIN  SOLUTION. 

RECALCULATE  ELEMENT  LOADS 

THE  REMAINING  FLAGS  ARE  NOT  USED  BY  PAFEC  LEVEL 
8.1. 

UNIFIED  CONSTITUnVE  MODEL  FLAGS  DESCRIPTION. 

30 

SUB-ITERATION  EXCEEDENCE  FLAG. 

31 

AUTOMATIC  (TRUE)  OR  MANUAL  (FALSE)  FLAG. 

32 

OUTPUT  PRINT  CONTROL  FLAG. 

33 

TABLE  PRINT  CONTROL  FLAG. 

34 

UNIFIED  CONSTITUnVE  MAIN  CONTROL  FLAG. 

35 

TURNING  POINT  INDICATOR  FLAG. 

36 

LOAD  INCREMENT  STORAGE  FLAG. 

37 

ITERATION  PRINT  FLAG. 

38 

DISPLACEMENT  PRINT  FLAG. 

39 

GAUSS  POINT  PRINT  FLAG. 

40 

STRESS  PRINT  FLAG. 

41 

CUMULATIVE  REACTION  PRINT  FLAG. 

42 

MAXIMUM  ITERATION  NON  CONVERGENCE  FLAG. 

43 

PREVIOUS  INCREMENT  TURNING  POINT  FLAG. 

44 

ELASTIC  UNLOAD  FLAG. 

45 ->>50 

PRESENTLY  NOT  USED 

/ 
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B.3:  Information  on  Non-linear  Data  Sets 

The  following  information  provides  a  description  of  the  variables  held  in  the 
linear  data  set  'NLSET'. 


POSITION 

CODE 

DESCRIPTION 

1 

* 

REFERENCE  STRESS  VALUES. 

2 

* 

STRESSES. 

3 

* 

REFERENCE  TOTAL  STRAINS. 

4 

* 

TOTAL  STRAINS. 

5 

* 

REFERENCE  ELASTIC  STRAINS. 

6 

* 

ELASTIC  STRAINS. 

7 

* 

REFERENCE  PLASTIC  STRAINS. 

8 

* 

PLASTIC  STRAINS. 

9 

O 

REFERENCE  CREEP  STRAINS. 

10 

O 

CREEP  STRAINS. 

11 

* 

REFERENCE  PLASTIC  STATE  VARIABLES. 

12 

* 

REFERENCE  PLASTIC  STATE  VARIABLE 

INCREMENTS. 

13 

* 

EQUIVALENT  STRESSES  AND  STATE  INDICATORS. 

14 

o 

LOAD  HISTORY  FLAGS. 

15 

O 

ROTATION  MATRIX  AND  LEFT  STRETCH  TENSOR. 

16 

i 

PRESENTLY  NOT  USED 

18 

19 

* 

REFERENCE  JACOBIAN. 

20 

* 

JACOBIAN. 

21 

* 

REFERENCE  TIME  STEP  DATA. 

22 

* 

TIME  STEP  DATA. 

23 

* 

REFERENCE  PLASTIC  BACK  STRESS. 

24 

* 

PLASTIC  BACK  STRESS. 

25 

* 

REFERENCE  SCALAR  VARIABLES. 

26 

* 

SCALAR  VARIABLES. 

27 

1 

PRESENTLY  NOT  USED 

50 

Code 

*  Data  sets  used  in  the  unified  constitutive  model  code. 
O  Data  sets  NOT  used  in  current  analysis. 


non- 
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B.4:  Created  Scratch  Module  Information. 

Several  scratch  modules  are  created  at  initialisation  time  to  hold  single  and  double 
precision  data.  The  identifier  numbers  are  stored  in  common  block  NEWFLOW  and 
witiiin  scratch  module  JPSRl  itself.  The  common  block  NEWFLOW  has  been  declared 
as: 

COMMON  /NEWFLOW/  JPSR1JPSR2,JPSR3,JPSR4JPSR5,JPSR6(10,2) 

A  description  of  the  contents  of  each  scratch  module  is  provided  for  reference  below: 


MODULE /PSRl  (Created  JPIOOO) 

This  module  holds  global  information  required  for  the  implementation  of  tire  unified 
constitutive  model  code.  A  description  is  provided  below: 


POSITION 

VALUE 

DESCRIPTION 

1 

0.0 

Total  time  at  any  point  in  solution. 

2 

UNIT 

Initial  delta  time  step.  Later  the  current  global  delta  time. 

3 

TOLTM 

Termination  time  of  solution. 

4 

0.0 

Copy  of  location  2  Used  when  sub-iteration  convergence  is  not  achieved. 

5 

IBLOK 

Block  size  to  increase  allowed  load  increments. 

6 

MODTI 

Module  number  holding  variables  for  Load-Time  interpolation  for  each  load 
case,  (Setup  JPlOlO,  used  IN  JPllOl,  JP1102) 

7 

NSEG 

Global  time  segment  pointer.  (Used  JPlOlO) 

8 

MODEL 

Module  number  holding  elastic  unload  time  information.  0  then  no  elastic 
unloads  exist.  (Setup  JPlOlO,  used  JP1106) 

9 

UNUSED 

10 

CONVTT 

Number  of  iterations  to  converge.  (Setup  NLllOO,  used  JPllOO) 

11 

TURTIM 

Storage  of  next  turning  point  time.  (Setup  JPllOl,  used  JP3920) 

12 

ABORFAC 

Abort  reduction  factor.  (USED  JPllOO) 

13 

AMAXTIM 

Maximum  time  step  allowed  (used  JPllOO) 

14 

TURDELT 

Delta  time  prior  to  a  turning  point.  Used  to  recover  delta  time  in  JPllOO. 

15 

ASTOR 

Delta  storage  time  value.  Controls  when  an  increment  is  written  to  disk. 

16 

CTIME 

Current  print  total  time,  (used  JP1103) 

17 

UNUSED 

18 

UNUSED 

19 

UNUSED 

20 

UNUSED 

MODULE  JPSRl  (Created  JPlOlO,  Used  JP1101,JP1105) 

This  module  holds  percentage  load  versus  time  values  for  each  load  case.  This  module 
is  created  from  the  expanded  incremental  module  and  does  not  exist  in  manual 
control.  The  created  module  is  printed  in  the  phase  9  output  file  so  that  the  load-time 
data  may  be  checked  by  tire  user. 
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MODULE  7PSR3 

This  module  is  currently  not  used. 


DSTO-TR-0529 


MODULE  JPSRi  (Created  and  Used  JP1120) 

This  module  is  used  to  store  die  previous  increments  cumulative  displacements.  This 
module  is  recovered  when  sub-iteration  convergence  has  not  been  achieved. 


MODULE  JPSR5  (Created  JPIOOO,  Used  JP1105,  JP3000,  p»3002) 

This  modules  holds  the  double  precision  control  data  required  by  kmodel.  This  module 
is  a  double  precision  copy  of  most  items  in  the  UNIFIED.CONTROL  module  231. 


MODULE  JPSR6{10,2)  (Created  JPIOOO,  Used  JP3000) 

This  array  holds  die  double  precision  material  data  module  number  in  location  1  and 
the  material  number  in  location  2.  The  module  number  referenced  is  a  double  precision 
copy  of  the  CONSTANTS.NEW.FLOW.LAW  module  228  less  the  material  number  for 
each  material  number.  Each  module  specifically  holds  aU  the  state  variables  associated 
with  each  temperature  for  a  particular  material  number.  A  maximum  of  10  different 
material  numbers  are  catered  for. 


MODULE  MODTI  (Created  JPlOlO,  Used  JPllOl,  JP1102) 

Module  number  holding  variables  for  Load-Time  interpolation  for  each  applied  load 
case.  Each  line  of  the  module  contains  the  items  described  below: 


ISEG 

TURTIM 

ITURN 

ACTIVE 

ALOAD 

4  Variables 


segment  in  load-time  curve. 

Next  turning  point  time  for  load  curve. 

Flag  to  indicate  that  this  curve  is  at  a  turning  point. 

Flag  to  indicate  that  this  load  case  is  active. 

Actual  extrapolated  load  value  at  the  current  time. 

Used  to  recover  load  and  accumulated  load  when  non-convergence 
recovery  is  required. 


MODULE  MODEL  (Created  JPlOlO,  Used  JP1106) 

Module  number  holding  elastic  unload  time  information  generated  by  negative  values 
of  time  in  the  TOTAL.TIME.VERSUS.FORCE.DATA  module  230.  If  no  elastic  unloads 
exist  then  this  module  number  is  initialised  to  zero  if  it  does  not  exist. 
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